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Cotunneling current through a resonant level coupled to either normal and superconducting or 
to two superconducting leads is studied for the domain of bias voltages, V , exceeding the supercon- 
ducting gap, 2A. Due to the on-site repulsion in the resonant level, cotunneling of an electron is 
accompanied by creation of a quasiparticle in a superconducting lead. Energy conservation imposes 
a threshold for this inelastic transport channel: Vc = 3A for N-S case and Vc = 4A for the S-S case. 
We demonstrate that the behavior of current near the respective thresholds is nonanalytic, namely, 

57*"(T/) oc (1/ - Vcf^'^ e{V~ Vc) and 5J™(y) <x {V - K)e (v - k) . Stronger anomaly for the 

S-S leads is the consequence of the enhanced density of states at the edges of the gap. In addition, 
the enhanced density of states makes the threshold anomalies for two-electron cotunneling processes 
in the Coulomb-blockaded regions more pronounced than for the N-N leads. 

PACS numbers: 73.23.-b, 73.23.Hk, 73.63.Kv, 74.50.-|-r 
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I. INTRODUCTION 

Early single-electron-transport devices^ were based on 
conducting grains containing large gate-controlled num- 
ber of electrons. A grain was coupled by tunnel barri- 
ers to two macroscopic leads. With number of electrons 
on the grain being large, superconductivity could be in- 
duced in the grain^'^ii^i^i^ii'^i'Siii^ upon lowering the tem- 
perature, while the leads remained cither normal*'^'^ or 
also turned into a superconducting state^ii^ii^iiiSiiiS. The 
focus of the early studies was the interplay between the 
two low-energji^i^i^iiiSiS scales, namely, the charging en- 
ergy and the superconducting gap. This interplay man- 
ifested itself in the Coulomb-blockade oscillations. On 
the theoretical side, different regimes of transport via 
superconducting grai n^^i^^'^^i^^ were studied for experi- 
mentally relevant situation of a grain containing many 
electrons. 

In the later experiments the grains have been 
replaced by much smaller few-electron quan- 
tum dots, based either on InA a^^^^'^^^^ or carbon 
nanotubea^ '^^i-^^'-^''i-^^'-^^i^°'^^'^^ . In these devices, there 
is no superconducting pairing of electrons on the dot. 
Rather, either one^^ or botbiii^^ii^ii^^ia^^^ 
leads are made of superconducting material. 

Interesting physics in the S-N-S junctions with super- 
conducting leads is due to the fact^ that the Andreev 
procesi 
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in these junctions gives rise to a rich subgap 
structure in the current-voltage characteristics^Si^i^. 
When the N-region is a small quantum dot (or a sin- 
gle resonant level) coupled by tunneling to the leads, this 
subgap structure is more pronounced22iMi2^. In addition, 
in the latter case the on-site interaction of two electrons, 
which in a small dot assumes the role of charging energy, 
becomes important^^^^^^^^^^. 

What makes the S-N-S structures with a resonant level 
as a N- region particularly interesting, is a delicate inter- 
play of a new energy scale, Kondo temperature, which is 
much smaller than the charging energy, and the super- 



conducting gap. This interplay is the focus of the very 
recent experimental studie a^^'^'^i^^'^^i^^'^^ . The results re- 
ported in Refs. il2i il3ij2l] suggest that subgap anomalies 
in differential conductance, G{V), at biases V — ±A, 
where 2A is the superconducting gap, are enhanced in 
the Kondo regime. Another intriguing observation made 
in Refs. IT^Il3ll2ll is that Kondo resonance leads to smear- 
ing of the conventional anomalies in G{V) a,t V = ±2A. 

Therefore, both theoretical^^^i^^^^^^^ and 

experimenta l^^'^'^1^^'^^1^^'^^ studies suggest that on-site 
repulsion affects the subgap structure in the conductance. 
However, it is commonly believed that for V > 2A there 
is no qualitative difference between the cases when su- 
perconducting leads are separated by a barrier or both 
coupled to a quantum dot. In the present paper we 
demonstrate that on-site repulsion manifests itself even 
for V > 2 A, leading to supergap anomalies in G{V). 
The underlying reason is that, at finite repulsion, in- 
elatic electron transitions between normal leads become 
possible's.. These transitions are accompanied by a quasi- 
particle excitations in the leads. When one of the leads 
is superconducting, the minimal energy of the excitation 
is 2A. Then, in order for electron tunneling from the 
normal lead to create the excitation in superconduct- 
ing lead, the bias should exceed Vc — 3A. Threshold 
for inelastic tunneling results in a supergap singularity, 
SG{V) oc {V — VcY^^ , in the N-S conductance, as demon- 
strated in Sect. Ill A. For the same reason, inelastic tun- 
neling between two superconducting leads has a threshold 
at Vc = 4A. We show that the supergap anomaly in the 
S-S transport has a step-like form, SG{V) oc Q{V — Vc) 
(Sect. IV), i.e., it is stronger than in the N-S case. This 
is due to the enhancement of the density of states at 
the edges of the superconducting gap. In addition, finite 
temperature, T, affects the S-S supergap anomaly only 
via the temperature dependence of A, whereas the N-S 
supergap anomaly is a universal function of {V — Vc)/T 
(Sect. Ill A.). In Sect. Ill B we also demonstrate that 
the enhancement of the density of states causes a sharp- 
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FIG. 1: (Color online) Inelastic correction to the lifetime of 
the localized state due to excitation of a quasiparticle across 
the gap is illustrated schematically. 



ening of the large-bias transport anomalies^ that involve 
two-electron transitions. 



II. ANOMALY IN THE LIFETIME OF A 
LOCALIZED STATE 

In order to illustrate how the on-site repulsion, U, gives 
rise to the anomalies in the conductance, G{V), we start 
from an auxiliary problem of the escape of an electron 
from the occupied localized state (LS) into a supercon- 
ductor. This situation is illustrated in Fig. 1. If the 
energy of the LS, Cd, lies above the upper boundary of 
the gap, 6(1 > A, then the population of the LS, which is 
occupied at time t = 0, decays with t as 



i{t) = exp(-rt). 



(1) 



where the decay rate, F, is given by the golden-rule ex- 
pression 



r(erf) = 7r7^t'o.9(ed). 
Here, 7 is the tunnel matrix element and 

z/og(e) = 1^0 



x/e2 - A2 



(2) 



(3) 



is the density of states in the superconductor. Eq. ^ 
applies when F is much smaller than A. Our main point 
is that, for large enough e^, there exists another inelas- 
tic channel of the electron escape into the continuum. 
Namely, the escape can be accompanied by excitation of 
a quasiparticle across the gap. This process leads to the 
threshold anomaly in the dependence F(ec;). The position 
of the threshold, e|j^'' , can be found from the following two 
conditions on the energy, e, of electron leaving the LS 



e>A, (ed-e)>2A. 



(4) 



The first condition ensures that the state into which elec- 
tron escapes is empty, while the meaning of the second 



condition is that the energy loss suffered by escaping elec- 
tron is sufficient to create a quasiparticle. From Eq. ^ 
we find the minimal value of ed 



3A. 



(5) 



Inelastic process is enabled by a finite U . To see this, we 
notice that there are two contributions to the amplitude 
of the process: 

1. An electron from the LS tunnels into the state e > A 
(i); another electron from the occupied state, e_, enters 
the LS (ii) , and subsequently tunnels into the empty state 
e+. These steps are illustrated in Fig. 1. 

2. Initial and final states are the same as in 1, while 
the intermediate steps (i) and (ii) are interchanged. As a 
result, after the first step, the LS is doubly occupied. In 
the absence of the on-site repulsion, the two amplitudes, 
1 and 2, would cancel each other identically. At finite 
[/, this cancellation does not happen. Note that, for 
large U ^ ed, the energy denominator corresponding to 
e_ — > ed contains U, so that the second amplitude can 
be neglected. 

The above reasoning is quite similar to that in Ref. (isl . 
where another inelastic process, occupation of the LS in 
the course of cotunncling between normal leads, has been 
considered. 

The amplitude, Al'^li^_ , of the three-step process in 
Fig. 1 is cx 7^^. Taking into account that the energies 
of the intermediate states are e and e -I- — e_ , the ana- 
lytical expression for this amplitude reads 
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Note, that this expression is valid when the states, e and 
£+ correspond to the opposite spin projections^, so that 
these states are distinguishable. On the contrary, for 
parallel spins of the states e and e+ the amplitude Eq. ([6|) 
vanishes'**' . 

The expression for inelastic correction to the rate, F, 
follows from Eq. ([6]) 



C30 00 — 

6r{ed) = y" dei^{e) J de+v{e+) j de_u{e_) 



X mtJ'S ed + e^-{e + e+) 



(7) 



It is seen from Eq. ([7]) that the argument of the 5- 
function turns to zero for — 3A at e_ = —A, and 
e = e+ = A. To establish the form of the anomaly near 



fid 



3A, we introduce the new variables 



A, 



(8) 



in Eq. ([7]). Now it is sufficient to set e_ = —A, and 
e = e-|_ = A in the denominator of Eq. ([7]), and re- 
place vie), v{e+), and i/(e_) by v^^fKplE, vq^/KJ^E^, 
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and V()yJ IS./ 2E- , respectively. Upon this replacement, 
Eq. ([7]) simplifies to 




X 5 ed-ty -{E + E+ + E^) 



(9) 



1/2 



the 



The above integral is proportional to ( 

numerical factor can be easily expressed through the sur- 
face area of the unit sphere. The final form of the thresh- 
old anomaly is the following 
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Ed 



1/2 



e 



Ac) 



(10) 



In deriving Eq. (jlOp we assumed that the intrinsic width, 
r, is much smaller than A. This guarantees that the 
relative correction 5T /T is small. The anomaly Eq. pO|) 
is much stronger than the threshold anomaly for two- 
electron ionization of the LS in Ref. [45|. The origin of 
this enhancement is the divergence of the density of states 
Eq. ^ at edges of the gap. 

In the above calculation we treated the states e_, e+, 
and e as electron states in a normal metal, and took su- 
perconductivity into account only via the energy depen- 
dence of the density of states, v{e). This is justified when 
the tunneling amplitude is calculated to the lowest order 
in the matrix element, 7. However the anomaly Eq. (jlOp 
emerges in the third order in 7. The proof of the validity 
of Eq. ^ for 5T{ed), starting from the BCS Hamiltonian, 
is presented in the Appendix. 




FIG. 2; (Color online) Origin of the anomaly at 14 = 3A in 
cotunneling between the N and S leads is illustrated schemat- 
ically. 



The elastic cotunneling conductance is given by 
4e2 TrVl 



Gel 



NS 



ttH (ed-F/2)2' 



(12) 



where we assumed (e^ — V/2) <^ V. The widths Tl.r — 
t^'^l.rIl r defined in a usual way; due to the en- 
hancement of the density of states in the superconduc- 
tor the width, Ffj, which enters into Eq. (fT2|l . becomes 

In order to calculate the inelastic correction, (5G™(1^), 
to the conductance, one cannot simply modify F^j ac- 
cording to Eq. (fTO| . This is because, in the course of 
cotunneling, the electron occupies the LS only virtually. 
The correct procedure of finding SG™{V) requires cal- 
culation of inelastic correction, SP"{V), to the current, 
taking into account that electron, transferred from the 
normal into superconducting lead, can excite a quasipar- 
ticle in this lead. Then we have 



III. SUPERGAP ANOMALIES IN THE N-S 
COTUNNELING 

A. Single-electron transport 

Passage of current from a metal to a superconductor by 
single-electron transitions, involving the LS, is illustrated 
in Fig. 2. Position, Vc, of the anomaly, at which the co- 
tunneling from the normal lead can be accompanied by 
creation of a quasiparticle in the superconducting lead, 
can be found from the similar reasoning as in Sect. I. The 
only difference is that electron enters the superconduct- 
ing lead with energy close the Fermi energy of the normal 
lead, so that 



K = 3A. 



(11) 



The magnitude of the anomaly is, however, weaker than 
for the electron escape considered in Sect. I. This is due 
to the fact that, while the energy of the LS is fixed to e^, 
the energy of the electron in the normal lead is simply 
restricted to the domain below V/2 - the Fermi level in 
the normal lead. 



Srr^{V) = ^ '^L'^liJ de,f{e, - V/2)l de g{e + V/2) 

-00 A-V/2 
00 -A-V/2 

X / de+g{e+ + V/2) f de_g{e. + V/2) 



A-y/2 



X \Al'X?5 ei+e_-(e + 6+) 



(13) 



where /(e) is the Fermi function, 
the transition amplitude (ei,e_) - 
Eq. ([6]) by an extra 7^, namely 



iLlR 



The expression for 
(e, e_|_) differs from 



(14) 



(fd - ei)(e - ei)(ed - e+) 

_^ 7L7| 

(e<i-e-)(e+-e-)(e<i-e)' 

As in Sect. I, in Eq. we had excluded the virtual 

states with doubly occupied LS. Two terms in Eq. ([M]) 
account for two different sequences in which the transi- 
tion (ei,e_) — > (e, e+) takes place. The first term corre- 
sponds to electron from the normal lead entering the LS 



FIG. 3: (Color online) The shape of the peak in the deriva- 
tive, dG/dV, of the N-S differential conductance is plotted 
from Eq. (|18|) versus dimensionless deviation x — {V — Vc)/T 



at the first step. The second term describes virtual oc- 
cupation of the LS by electron from superconductor with 
energy e_ at the first step, followed by its escape into e+ 
and subsequent cotunneling of electron from the normal 
lead. Note, that there is no analog of the second contri- 
bution in the amplitude Eq. This is because Eq. ^ 
describes the process in which the LS was occupied in 
the initial state. 

In order to extract the anomaly, upon substituting 
Eq. pi)) into Eq. we introduce the new variables 

Ei=ei^V/2, E = e + V/2-/^, (15) 
E+ = e+ + y/2 - A, E^ = -e_ - V/2 ~ A. 

For bias, V, close to Vc — 3 A, characteristic values of 
E, El, E+, and E- are much smaller than A. This al- 
lows to set ei = Vc/2 and e — ^ ~Vc/2 -I- A in the 
denominators of Eq. p4|) . We can also use the near-gap- 
edge asymptotes for the densities of states in the super- 
conducting leads. After these simplifications, Eq. (E 
assumes the form 



2^/2 g 



7T^ h 

TlTI (AT)3/2 



(16) 



[(e,~3A/2)((ed + 3A/2)2-A2)]- 



V~Vc 
T 



where the dimensionless function F of a single argument, 
iV~'Vc)/T, is defined as 



V^Vr. 



T 



dEif{Ei) fdE fdE. 
T3/2 



Ve 





X S V~Vc + Ei- (E + E^ 




Note, that the three- fold integration over E, E+, and 
E- has already been carried out in Sect. 1. It yields 




FIG. 4: (Color online) Stability diagram for transport be- 
tween N and S leads via a localized state. White region cor- 
responds to the sequential tunneling transport. Horizontal 
blue lines, V = ±3A, correspond to the single-electron super- 
gap anomaly, illustrated in Fig. 2. Red lines, V — 2/3(ed± A) 
and V — 2/3(ed ± 2A), are the positions of the two-electron 
resonance. Subgap resonances at \V\ < A lie in the shaded 
region. 



27r{V~Vc-Ei)^^^e 



V-Vc + Ei 



As a result, the bias 



dependence of J/'" is given by a single integral 

oo 

'V-Vr' 



T 



2n / dx 



exp [a 



1 



(18) 



Inelastic correction , SG™ {V), to the differential conduc- 
tance is thus described by the derivative, dF/dV. The 
asymptotic behavior of (5G'"(y) at low, {V - Vc) < T, 
and high, (V—Vc) ^ T, temperatures can be easily found 
from Eq. ([T5)) . We present the results for a dimensionless 
ratio, SG"^ /G'^'', of inelastic and elastic contributions to 
the conductance 



~G^ 



r FV^' ("l^-^c)»T, 

37r[(£rf+3A/2)2-A2j V ^ ^ 

, r pV?, {V-Vc)«T, 



(19) 



where a = 2~^/2 dxx^/'^ / cosYv' x ~ 0.536. It is seen 
from Eq. (fT9|) that, a,t V > Vc, differential conductance 
acquires a correction cx (V—Vc)^^'^. Correspondingly, the 
second derivative, (Pl/dV"^, has an asymmetric peak of a 
width '-^ T centered aiV = Vc- The shape of the peak is 
given by the second derivative of the function F. In Fig. 3 
this derivative, calculated numerically from Eq. (|18p . is 
plotted versus dimensionless deviation, {V — Vc)/T. 
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B. Two-electron transport 

1. Ionization of the LS 

In terms of the Coulomb blockade stability diagram in 
the (cd, V) plane, Fig. 4, the anomahes at V = ±3 A 
correspond to horizontal lines, which start from the 
points (-A/2, 3A), (3A/2, 3A), and (-3A/2, -3A), 
(A/2, — 3A). These lines extend into the blockaded re- 
gion. In Ref. it was demonstrated that, without su- 
perconductivity, there exists an additional weak struc- 
ture within the Coulomb blockade diamond, along the 
lines V — ±26^/3. The origin of this structure is the 
two-electron ionization of the LS, namely, the process, in 
which one electron from the left lead is transferred to the 
right lead while the other electron from the left lead occu- 
pies the LS. The position of the boundary, V = 2e(i/3, ex- 
presses the threshold for this two-electron transfer, which 
follows from the energy conservation. In this subsection 
we point out that, in the presence of the superconductiv- 
ity, the boundaries for two-electron ionization are modi- 
fied in an asymmetric fashion. For positive bias, > 0, 
the boundaries are located at 

T/+(6,)=±^(^6,-^)+A, (20) 

while for negative bias they are located at 

F_(e,) = T^(e<i + y)-A. (21) 

These modified boundaries are shown in Fig. 4. More 
importantly, as we demonstrate below, superconductivity 
leads to the strengthening of the ionization anomaly. The 
underlying mechanism for this strengthening is, again, 
the enhancement of the density of states at the bound- 
aries of the gap. 

Energy dependence of the density of states can be eas- 
ily incorporated into the expression from Ref. [4^ for ion- 
ization rate. Consider first the situation when the initial 
states of two electrons with energies ei and £2 are in the 
normal lead, while one of the finite states (with energy 
e) is in the superconducting lead and the other is on the 
LS. The ionization rate for T = is given by 



V/2 V/2 00 

dei J de2 J de g{e + V/2) 

-00 A-V/2 



1 



(ed-ei)2(ed-e2)2 



ei - £2 



(22) 



Near the threshold, V = V+{ed), one can set ei = £2 = 
V+/2 in the denominator of Eq. ((22|) . Upon measuring 
the energies £1, £2, and e from their respective bound- 
aries, as in Eq. (llSp . we can simplify Eq. (|22p to 



r?r«A2 



(2^)2(e,-y+/2)' 



V-V+ 



(23) 




FIG. 5; (Color online) Shapes of the anomalies in cfl/dV^ 
near V = V+ and V — V- versus dimensionless deviations 
X = iy — V±)/A calculated, respectively, from Eqs. (|25[) 
(black line), and (|3Up (red line). Both curves approach the 
value 9/4 at x 00. For small x, d^H+/dx'^ diverges as 
33/2/4^1/2^ while the a; = value of d^H-fdx^ is 277r/32. 
Note that H+{x) and H-(x) are zero for 2; < 0. 



where the one-parameter function, H+, is defined as 

00 00 

TT / ^ [dEi [dE2 /3 , , \ 



xe 



3 El E2 

2""" A"" a: 



(24) 



This integral is easily calculable. Its analytic form is 



H+(:e) 



3x + 2 
8 

-Im 
2 



(25) 



3 / 9 

-X + 1 + \ -x'^ + 3x 

2 V 4 



The large- x and small- x asymptotes of H+ are 
\/3x3/2, x<l. 



H+(x) = 



9 ^2 



(26) 



x"^, a; > 1, 



Consider now V < 0. The ionization rate is given by 
the expression 

-V/2-A 



-V/2-A 



de2g{e2 + V/2) J de 

v/2 



i 



£d + £ - £1 - £2 



(£<i-£l)2(£d-£2)2 



(27) 

which differs from Eq. by additional density of su- 
perconducting states in the integrand. When the bias 
voltage is near the critical, V = V^, one can replace the 
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values of ei and £2 by their boundary value —V-jl — A, 
in the denominator of Eq. (|27p . Then the ionization rate, 



^ion^ can be expressed as 



'[V) 



(2^)2(e, + F_/2 + A)4 



H 



V- - v 



(28) 



where we have absorbed all the integrals of Eq. ([27|) into 
the new one-parameter function, H_, defined as follows: 

H- (-) 9iE, + A)|^ giE, + A) 


El 



xe 



2"^ 



E2 
a" 



(29) 



One of the integrations in Eq. can be performed 
explicitly. The final form of the function H„ is the fol- 
lowing 



3x/2 



[z + l)v/3x/2 - z ^ixjl + 2- z 
+ 2z 



(30) 



The easiest way to find the behavior of H_(x) at large 
and small x is to set, respectively, g{t) = 1 and 5(e) = 
■y/A/2e in the integrand of the definition Eq. ([29]) . This 
yields 

fx (1 + f^x) , X « 1, 
H_(x) = <( ^ ^ (31) 

^ X , X 1 , 

Upon populating the LS, the electron rapidly, within the 
time {Tl -\- ^r)~'^, escapes either to the left or to the 
right lead. In terms of contributions to inelastic cur- 
rent these two channels of escape are differentia. For 
escape to the left, the net charge transfer is e, while 
for escape to the right, it is 2e. As a result, the in- 
elastic contribution to the current is equal to 5P'^{V) = 
2er„„ {V){2Tr + Tl)/{Tl + Yr)- The threshold behavior 
of 51"'' near V+ and F„ is determined by the functions 
i7+ and respectively. As seen from Eqs. (f26| and 
([31]) , these behaviors coincide when V — V+, V ~ are 
much bigger than A. This is natural since for large de- 
viations from the thresholds, superconducting gap drops 
out from 5J'". Note, however, that in the immediate 
vicinities of V+ and V- , the threshold behaviors are dif- 
ferent, namely, SP"^ is more singular near V- than near 
V-^. The origin of this asymmetry is that the inelastic 
process N S involves only one state near the supercon- 
ducting gap, while the inelastic process S N involves 
two such states. Without superconductivity, threshold 
anomaly at V = zL2ed/3 shows up in the third derivative 
of the current with respect to V. Our results, Eqs. ([22]) 
and ([28]) . suggest that, within the interval ~ A from 
the thresholds V+, the singularities of current are 



more pronounced: they show up already in the second 
derivative (Pl/dV^. This is illustrated in Fig. 5, where 
d'^U+fdV^ and d'^B.^/dV^ are plotted. 

In fact, the singular behavior of inelastic current near 
V — V- shows up already on the level of differential 
conductance, JG" = dSP'^/dV, as a step oc e(VL - F). 
Combining Eq. (|28|) and Eq. (j31|) . we get the following 
magnitude of the step 



3e^ 



16nh {ed + V-/2 + A) 



-e{V.-V). (32) 



2. Two-electron tunneling 

As seen from Eqs. and (P^ . the relative correc- 
tion, SP"{V)/P'', to the elastic current due to ionization 
of the LS, changes on the scale ~ A; the magnitude 
of correction a,t V ^ A being ^ F^A/e^ ^ 1. Although 
small, this correction is distinguishable by virtue of its 
threshold dependence on bias. Indeed, both H_|_ and H_ 
are zero for V < V+ and V > V-, respectively. Another 
fact that distinguishes the transport at biases near V+ 
and V- is that the inelastic current, SP"-{V), has a pre- 
cursor with singular dependence on deviation V —V± and 
on the temperature T. The origin of this precursor— is 
direct cotunneling of two electrons via the LS. This pro- 
cess differs from ionization of the LS, since, in course of 
this two-electron cotunneling, the LS is populated only 
virtually. As a result, the corresponding contribution to 
the current, I^{V), contains extra power F^ (or Tn). 
On the other hand, this contribution is more singular in 
deviation, V —V±, and has a peculiar T-dependence. As 
all other corrections to the elastic cotunneling calculated 
above, I^{V) is enabled by a finite on-site repulsion. The 
golden-rule expression for l'^{V) 



ft(2^)2(erf-y+/2)2(14 + A)2 

00 

X / dei /(ei - V/2) [ de2 /(es - V/2) (33) 



-00 
00 



— 00 
00 



dEi g{Ei + V/2) / dE2 



-V/2+A 



+ £2 ~ El - E-, 
(ed - E^r 



contains energy denominators that correspond to virtual 
states; in these states the LS is occupied by first and 
then by second tunneling electron. Note, that the 5- 
function in Eq. ([33)) ensures conservation of the total 
energy, ei -\- €2, of two electrons in the initial and final 
states, while individual energies get redistributed. Sen- 
sitivity of 5l'^{V) to V = V+ comes from the domain 
of integration in Eq. ((33)l with Ei near the Fermi edge. 
El w ~V-^/2 + A, and E2 ~ e^- For this reason, the 
nonresonant energy denominators are extracted from the 
integrand of Eq. 
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As before, in Eq. p8|) we used the near-gap-edge asymp- 
tote of g{e), so that Eq. ([38|) apphes in the interval 

rL,rfl« |i/-y+|,r« A. 

The four-fold integration in Eq. (l34|) can be reduces to 
a single integral by using the Fourier representation for 
the 6- function and the fact that the Fourier transform 
of the Fermi function is equal to /{lo) = 7r/sinh(7rti;). 
We will present the result for the second derivative, 
/ dx'^ , which describes the near-threshold behavior 
It reads 



of (fll^/dV 



FIG. 6: (Color online) Shapes of the anomalies in Sl'± /dV'^ 
near V = V+ and V = V- versus dimensionless deviations 
X = {V ~ V±)/T calculated, respectively, from Eqs. p4)) 
(black line), and (I39|l (red line). 



In order to capture the dependence of 61^ on {V — 
V+) and T, we introduce the dimensionless function, K+, 
defined as 



dx^ 



' 25/2 



(36) 



ds- 



«5/2 



sinh (tts) 



cos 



sm I — sx 



K+{x) 



dzi 



dz-). 



e^i +1 J + 1 7 ^ 



/ ^^£4 d'^K+fdx'^ at large negative a; is oc \x\ 



dz4 

72 



Consider first the limiting case of vanishing T. To re- 
alize that the temperature drops out from the expres- 
sion Eq. we notice that the asymptotic behavior of 



-3/2 



This yields 



xS 



-X + Zl + Z2 - Z3 - Z^ 



Then 5l']^{y^ can be presented in the form 



e TITr{2Tr + Tl) [AT 
h{27r)^{ed-V+/2y{V+ + A)^ \ 2 



T 



d'^61^'' /dV^ oc (T^ — V+)~^^^. The divergence is stronger 
than 1/{V — V+) in Ref. \^ Another remarkable fea- 

(34) ture of d^I^/dV^ is that, at finite T, it exhibits a fine 
structure. This is seen from Fig. 6, where the function 
d'^¥i+ /dx^ is plotted. Asymptotic behavior of (i^K+ /dx'^ 
at large positive a; is cx exp(— 3a;/2). This suggests that 

1/2 for V > Vj^ ionization current dominates over 51^ . 

Calculation of the two-electron current, 51'^ {V), near 
V = V- IS quite similar to Eqs. ([33]) and (l38|) Namely, 

(35) the golden-rule expression 



7lT\{Tr + 2rL) 



-y/2-A 



-y/2-A 



h {2nf{ed + V-/2 + ^f{V- + A)^ 



dE^ [1 - f{E^ - V/2)\ / dE2 



J dei g{ei + V/2) J de^ 5(^2 + V/2) 



d(^ei+e2 -Ei-E^^ 



[ed - E2Y 



(37) 



is cast into the form 



and final states 



6ll^{V) 



+ 2r^)A 



00 00 



/i 87r2(ed + F_/2 + A)2(F_ + A)2 
'V- - V 



K 



T 



dzi f dz2 



(38) 



—00 —OQ 



X 6 



-X + Zi + Z2 + Z3 + Z4 



(39) 



The dimensionless function K_(a;) is a four-fold integral 
over the electron energies (in the units of T) in the initial 



Three out of four integrations in Eq. (|39p . over zi, Z2, 
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FIG. 7: (Color online) One of the possible inelastic channels in 
the S-S transport is illustrated schematically. Electron trans- 
fer is accompanied by creation of the excitation in the left 
lead. 



and Z4, can be carried out explicitly. Then we get 

oc 

■ 9 = ^ / ds^— — - cos -sx . (40) 

dx^ 4 J sinh(7rs) V2 / 



Second derivative, cPK.-/dx^ , is plotted in Fig. 6. It 
shows that d^I'^/dV'^ also exhibits a fine structure at 
(V- -V)r^T. The low-T behavior of d'^P^dV'^ is even 
more singular than that of d^I^ /dV^. This follows from 



the large- cc asymptote, oc of the integral Eq. (|40p . 

Thus, the below-threshold behavior of d^I^ /dV^ is oc 

\/{V-~vY. 

IV. ANOMALY IN THE S-S COTUNNELING 

Energy diagram for transport between two supercon- 
ducting leads via an LS is shown in Fig. 7 for bias 
V > 2 A. Similarly to the case of normal and supercon- 
ducting leads, electron cotunneling can be accompanied 
by excitation of a quasiparticle across the gap. It is easy 
to see from Fig. 7 that the threshold bias for this process 
is 

Vc = 4A. (41) 

The difference from the N-S case is that, at the threshold, 
electron tunnels from the edge of the gap rather than 
from the Fermi level of the metal. A more significant 
qualitative difference from the N-S geometry is that a 
quasiparticle can be excited in both leads. Besides, as we 
will see below, the anomaly is stronger in the S-S than in 
the N-S case. This is due to the divergence of the density 
of states in both leads. More specifically, instead of the 
four-fold integral Eq. , the near-threshold expression 
for inelastic contribution to the current reads 



ft^3[(,^ + 3A)(e2_A2)] 



V/2-A 

2 I deM^i - V/2) I de g{e + V/2) / de+g{e+ + V/2) 



A-y/2 



A-V/2 



X /de_5(e_ + y/2)<5 ei +e_ - (e + e+) 



(42) 



Similarly to Eq. (fT5|) . in order to calculate the integral 
Eq. (|42| . we introduce the same variables E, 
as in Eq. and also Ex = -(ei - V/2 -I- A). Upon 

taking the near-gap asymptotes for the density of states, 
Eq. (H^ assumes the form 



OO ^ OO 

dEi fdE 



^7r3[(ed + 3A)(e2_ A 

OO 

dE+ fdE^ 



V y/Ei J VE 




^E+J JE. 

" 



V -Vc- Ei+E + E++E^ 



(43) 



After rescaling all variables to {V — Vc), this integral re- 
duces to the surface area of a unit sphere in four dimen- 



sions, and we obtain: 



5/™(X/) 



^127r2[(erf + 3A)(e2-A2)] 
x{V-Vc)Q \v-Vc 



(44) 



Contribution Eq. describes cotunneling accompa- 

nied by excitation of a quasiparticle in the right lead. 
Similar calculation for inelastic channel, with excitation 
of a quasiparticle, as depicted in Fig. 7, results in 



5/™(F) = 1 



TrTIA' 



^127r2[(e,-3A)(e2~A2)] 
x{V-Vc)e \v-Vc 



(45) 
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Here we would like to emphasize that both calculations 
leading to Eqs. and ([iS]) take into account that quasi- 
particle can be created at the first as well as at the last 
step of the cotunneling process, and corresponding am- 
plitudes interfere^ as in Eq. (jl4p . Taking this interfer- 



ence into account, results in the extra factor ~ ^^/e^ 



in Eqs. (|44l) and (|45p . Obviously, the threshold anomaly 
Eq. (I44p in the current results in the jump in the In- 
dependence of the differential conductance. Within a 
numerical factor and assuming ^ A, the magnitude 
of the jump can be presented as 



5G' 



5G' 



(46) 



Here the sum T\ 



r|j, accounts for the contributions of 
the two channels of inelastic current, mentioned above. 
We note, that the step Eq. ([i^ is abrupt; its temperature 
smearing is oc exp(— A/T) rather than ~ T, as in the case 
of tunneling between N and S leads. 

Overall, the stability diagram for superconducting 
leads differs from Fig. 4 in two respects. Firstly, the 
positions of the supergap anomalies are Vc = ±4A. 
Secondly, the stability diagram is symmetric with re- 
spect to V ^ —V. Namely, the boundaries of the two- 
electron ionization anomaly in this case are located at 
V* = ±|(ed-F3A). Regarding the "strength" of two- 
electron anomaly, the threshold behavior of the differen- 
tial conductance can be found from calculation similar to 
Eqs. (P^ . (P7)l . and within a prefactor yields 



SG'^'iV) cx {V~V*) 



-1/2 



(47) 



i.e., the threshold behavior is more singular than 
Eq. ([5^ . Again, the divergence of i5G™ is limited by 
iy —V*)^ ^L,R rather than by temperature. 



V. CONCLUDING REMARKS 

Let us list the assumptions adopted in the above con- 
sideration: 

(i) energy position, e^, of the LS is well outside the su- 
perconducting gap. A; 

(ii) on-site repulsion, C/, is the largest energy scale, U ^ 

(iii) the widths, F^, F^j, are the smallest energy scales, 
so that 



VlJr < a < < t/. 



(48) 



One of the consequences of Eq. (|48|) is that the Kondo 
temperature, Tk csc exp[— 7r|ed|/2(Fi -I- Fj^)], is much 
smaller than A. This means that the Kondo effect will 
not developed fully, but rather manifest itself as an en- 
hancement oc ln~^A/rR- of the conductance at small 
bias. 

Recent experimental papers Refs. [HIllIllIlUlllll 

are focused on the domain of parameters Tk ~ A, where 



the two prominent regimes of transport compete with 
each other. This competition is due to the fact that an- 
tiparrallel spins of electrons in the Cooper pairs cannot 
mediate the spin-flip processes that are responsible for 
the Kondo effect. Experimentally, in the case of normal 
leads, the Kondo effect manifests itself on the stability 
diagram in the (e^, V) plane as enhanced zero-bias con- 
ductance in the valley td < 0, where LS is occupied. It 
has no effect on the valley > 0. On the other hand, 
with superconducting leads, conductance is suppressed 
in the entire domain of biases V < 2A in both valleys. 
A non-trivial result of interplay between the Kondo ef- 
fect and superconductivity is that the peaks at V = ±A 
emerge in the Kondo valleys, whereas the conventional 
peaks a.t V = ±2A are suppressedi^J^i^. This implies 
that the Andreev transport process is facilitated by the 
Kondo resonance. Conversely, in the non-Kondo valleys, 
the peaks V = ±A do not show up, while V — ±2A- 
peaks are strong and exhibit a well-known threshold be- 
havior, reflecting the BSC density of states. 

In the present paper we predict additional anomalies 
both outside the Kondo regime and above the gap. Nev- 
ertheless, the origin of the new anomalies is intimately 
related to the Kondo physics. To clarify this relation, we 
recall that, in a bulk metal with magnetic impurities the 
energy exchange between electrons is possible even with- 
out direct electron-electron interaction. This was first 
demonstrated by Kaminski and Glazman in Ref. [13. Ob- 
viously, such an exchange is impossible in the case of 
non-magnetic impurities. The reason is that the mech- 
anism, which is responsible for an impurity (LS) being 
magnetic, is a finite on-site repulsion, U. As a result, the 
interaction between two electrons in metal, leading to the 
energy exchange, takes place when they virtually visit the 
LS. The energy exchange occurs between electrons with 
opposite spins, and in the case of magnetic im pur ity, in- 
volves spin- flips'*''. Thus the mechanism Ref. |47j repre- 
sents the most elementary manifestation of the Kondo 
physics, and even does not require the presence of the 
Fermi sea. 

As was demonstrated in Ref. |4^ the mechanism^ can 
be extended to the transport between two normal leads, 
coupled to the LS. Then, for two electrons tunneling be- 
tween the leads, the magnitude of the energy exchange is 
limited by the applied bias, V. This leads to the anoma- 
lies in conductance at y = ±2ed/3. The main message 
of the present paper is that, in the case when one or both 
leads are superconducting, the gap, 2A, sets the thresh- 
old for inelastic process of one-electron transfer accompa- 
nied by a quasiparticle excitation in the superconducting 
lead. The ensuing anomalies at Vc = ±3A (for N and S 
leads) and at Vc = ±4A (for S-S leads) are independent 
of the gate voltage, e^. The anomaly near V — V is not 
smeared by temperature and manifests itself as a sharp 
peak in the second derivative (Pl{V)/dV^. Although the 
papers on transport through Coulomb-blockaded dots re- 
port the data on first derivative, i.e., the differential con- 
ductance, dI(V)/dV, the second derivative was previ- 
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ously measured for single-electron transport through a 
molecule^. In Ref. the second derivative was re- 
quired to resolve a fine structure in the /(y)-dependence, 
related to the vibrational satellites. 

As a final remark, we note that higher-order, in pa- 
rameters, Fl/A, Ffl/A processes will lead to anomalies 
at even larger biases due to creation of more than one 
quasiparticle by a tunneling electron. For the case of the 
S-S leads, additional anomalies can be expected at biases 
Vc "'^ = 2A(2 -|- n). Estimate for the behavior of inelastic 

current {V — Vc^^) <^ A can be easily found by extend- 
ing the four-fold integral in Eq. (|43p to higher n. This 

yields: 5P^'{V) {V ~ yi"V+ie(y - V^). 
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for the upper branch in Eq. (|49p and 



V2 



1 - 



1/2 



(51) 



for the lower branch in Eq. (|49p . Then, when perform- 
ing summation over states corresponding to, say, upper 
branch, one has to take into account contributions oc it(^) 
and (X u(— ^), since they correspond to the same energy 
e = a/A^ -f Now the fact that the main contribution 
T{ed) Eq. ([2]) to the lifetime does not contain coherence 
factors can be formally interpreted as a consequence of 
the identity u^{i) + v?{-i) = 1. 

Turning to the third-order amplitude Eq. ([S]), the cor- 
rect way to write one particular contribution to A\'^X- 
is 



VI. APPENDIX 

The fact that the superconductivity manifests itself in 
the expression Eq. ^ for the tunneling rate T{ed) only 
through the density of states Eq. ([3]) is well known. How- 
ever, it is not obvious that higher-order, in the tunnel 
matrix element, 7, corrections to r(erf) can be expressed 
solely through (?(e), and do not contain coherence fac- 
tors. Indeed, in our calculations we treated the ampli- 
tude, A\"^^_ , Eq. ([S]) as a number determined only by 
the energies 



(49) 



of initial and final states in superconductor^ and ignored 
the fact that the real amplitude contains contributions of 
positive and negative bare energies ^. This contributions 
enter into the amplitude with different weights, namely 



1 

71 



1/2 



(50) 



e+) = 



(52) 



e.-VAM^ 



-./A2 



-e-x/A^+F 



Then the correction 5T{ed) is, actually, the sum of all 
possible contributions, i.e.. 



jF(e,)(x E E E 

5_<0 S<0 {+<0 



(53) 



From Eq. ([55]) it becomes apparent that coherence fac- 
tors in the numerators of eight contributions can be 
combined into the product [u^(^_) -I- u^(— ^_)] [u^(^) -|- 
M^(— ^)][m^(^+) + u^(— ^+)], which is an identical unity. 
Note, that this conclusion rests on the assumption that 
the matrix element, 7, is independent of ^. 
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